# -*- coding: utf-8 -*-
###################################################
#    This file is part of blockIMH.
#
#    blockIMH is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License, or
#    (at your option) any later version.
#
#    blockIMH is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with blockIMH.  If not, see <http://www.gnu.org/licenses/>.
###################################################
#! /usr/bin/env python
from __future__ import division
from numpy import random, exp, sqrt, power, log, newaxis, zeros, \
        zeros_like, float32, int32, repeat, minimum, ones, mean, \
        average, sum, var, array, matrix, transpose, dot, arange
from numpy.linalg import inv, det
import rpy2.robjects as robjects
from scipy.stats import norm


r = robjects.r
r("""
library(MASS)
data(Pima.te)
#a <- glm(Pima.te[,8] ~ Pima.te$age+Pima.te$glu+Pima.te$bp+Pima.te$ped+Pima.te$npreg+Pima.te$bmi+Pima.te$skin,family=binomial("probit"))
#a <- glm(Pima.te[,8] ~ 0 + Pima.te$glu + Pima.te$bp + Pima.te$ped, family=binomial("probit"))
#a <- glm(Pima.te[,8] ~ -1 + Pima.te$glu + Pima.te$bp, family=binomial("probit"))
a <- glm(Pima.te[,8] ~ -1 + Pima.te$glu + Pima.te$bp + Pima.te$ped, family=binomial("probit"))
covmatrix <- summary(a)$cov.unscaled
estim <- a$coeff
y <- Pima.te$type
X <- Pima.te[,c(2, 3, 6)]
#X <- Pima.te[,c(2, 3)]
#X <- Pima.te[,1:7]
""")

estimMLE = array(r["estim"])
parameterdimension = estimMLE.shape[0]
scalefactor = 3
covarianceMLE = scalefactor * array(r["covmatrix"])
precisionMLE = inv(covarianceMLE)
observations = array(r["y"]) - 1
nbobservations = observations.shape[0]

XT = array(r["X"])
X = transpose(XT)

#print dot(X[0,:], estimMLE[1:8])

priorcovariance = nbobservations * inv(dot(transpose(X), X))
priorprecision = inv(priorcovariance)

########################
def logprior(theta):
    return -0.5 * dot(transpose(theta), dot(priorprecision, theta))

def loglikelihood(theta):
    z = dot(X, theta)
    z = observations * log(norm.cdf(z)) + (1 - observations) * log(norm.cdf(-z))
    return sum(z)

########################
def rproposal(size):
    proposals = random.multivariate_normal(estimMLE, covarianceMLE, size = size)
    return proposals

def logdproposal(theta):
    return -0.5 * dot(transpose(theta - estimMLE), dot(precisionMLE, theta - estimMLE))

def fastomegasingle(theta):
    results = logprior(theta) + loglikelihood(theta) - logdproposal(theta)
    return results

def fastomega(theta):
    results = zeros(theta.shape[0])
    for i in xrange(theta.shape[0]):
        results[i] = fastomegasingle(theta[i, :])
    return results


#print fastomega(proposals)
#x_starting = estimMLE
#print x_starting
#print fastomegasingle(x_starting)
#print covarianceMLE
#raw_input()


#print "simple IMH"
#
#def simpleIMH(T, x_starting):
#    wcurrent = fastomegasingle(x_starting)
#    chain = zeros((T + 1, parameterdimension))
#    chain[0, :] = x_starting
#    currenttheta = x_starting.copy()
#    uniforms = random.uniform(size = T)
#    proposals = rproposal(size = T)
#    wproposals = fastomega(proposals)
#    nbaccepts = 0
#    for t in xrange(T):
#        if t % 10 == 0:
#            print "iteration %i" % t
#        wproposal = wproposals[t]
#        acceptratio = wproposal - wcurrent
#        u = uniforms[t]
#        acceptation = (log(u) < acceptratio)
#        nbaccepts += acceptation
#        if acceptation:
#            chain[t+1, :] = proposals[t]
#            wcurrent = wproposal
#            currenttheta = proposals[t]
#        else:
#            chain[t+1, :] = chain[t, :]
#    print "accept rate %.4f" % (nbaccepts / T)
#    return {"chain":chain, "proposals": proposals}
#
#T = 10000
#indices = arange(start = 0, stop = T, step = 1)
#imh = simpleIMH(T, x_starting)
#robjects.r(""" par(mfrow = c(%i, 2)) """ % parameterdimension)
#for i in range(parameterdimension):
#    x = imh["chain"][indices, i]
#    rx = robjects.FloatVector(x)
#    robjects.r.plot(rx, xlab = "", ylab = "", type ="l")
#    robjects.r.hist(rx, xlab = "", ylab = "", main = "", col = "red", prob = 1, nclass = 100)
#raw_input()
#
#
#
